
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <mpi.h>
#include <omp.h>
#include "jacobi.h"
#include "matrixparser.h"


double max2(double *vdif, double *vx, int size)
{
    double max,max2;
    /*int i;
    double temp,temp2;*/

    max = max2 = 0;

    /*for(i = 0; i < size; i++)
    {
        temp = fabs(vdif[i]);
        temp2 = fabs(vx[i]);
        if(temp > max)
            max = temp;
        if(temp2 > max2)
            max2 = temp2;
     }*/

    omp_set_num_threads(1);
    #pragma omp parallel
    #pragma omp sections
    {
        #pragma omp section
            max = diff(vdif,size);
        #pragma omp section
            max2 = diff(vx,size);
    }

    return max/max2;
}

double diff(double *diff, int size)
{
    int i;
    int max = 0;
    double temp;

    for(i=0;i<size;i++)
    {
        temp = fabs(diff[i]);

        if(temp>max)
            max = temp;
    }

    return max;
}

double *solve_mpi(int argc, char *argv[],double *a, double *b, double *x, int m_order, int row_test,int max_it, double eps)
{
    unsigned int m_order_pow_2 = m_order * m_order;
    int j = 0;
    int i = 0;
    int it;
    int lines_per_task;
    int lines_tail;
    int index;
    int indexAdjust;
    double globalMax;
    double *recvbuf;
    double *bbuf;
    double diagonal_temp;    
    double *x_res;
    int tam;
    int tail;
    double *xTmp = malloc(m_order * sizeof(double));
    int t;
    double max = 0.0;
    double x_max = 0.0;
  
    int loop = 1;

    /*mpi declarations*/
    int rank,source,numtasks,sendcount,recvcount;
    
    /*new*/
    double *dif = malloc(m_order*sizeof(double));
    double *dbuf;    

    MPI_Init(&argc,&argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &numtasks);

    if(numtasks > m_order)
    {
        lines_per_task = m_order;
        tam = m_order;
        lines_tail = 0;
    }else {
        lines_per_task = m_order/numtasks;
        lines_tail = m_order%numtasks;
        tam = lines_per_task * m_order;
    }

    tail = lines_tail * m_order;
    x_res = malloc(m_order * sizeof(double));
    recvbuf = malloc(tam * sizeof(double));
    bbuf = malloc(lines_per_task * sizeof(double));

    /*new*/
    dbuf = malloc(lines_per_task*sizeof(double));

    source = 0;
    sendcount = tam;
    recvcount = tam;

    /*Send lines to tasks to process*/
    MPI_Scatter(a,sendcount,MPI_DOUBLE,recvbuf,recvcount,
             MPI_DOUBLE,source,MPI_COMM_WORLD);
    MPI_Scatter(b,lines_per_task,MPI_DOUBLE,bbuf,lines_per_task,
             MPI_DOUBLE,source,MPI_COMM_WORLD);

    /*new*/
    MPI_Scatter(dif,lines_per_task,MPI_DOUBLE,dbuf,lines_per_task,
             MPI_DOUBLE,source,MPI_COMM_WORLD);

    /*Process lines divided for each process*/

    for(i = 0; i < lines_per_task; i++)
    {
        index = lines_per_task*rank;
        index = index +((m_order+1)*i);
        diagonal_temp = recvbuf[index];
        recvbuf[index] = 0;
        indexAdjust = i*m_order;

        for(j = 0; j < m_order;j++)
        {
            if(j != index)
                recvbuf[j+indexAdjust] /= diagonal_temp;
        }

            bbuf[i] /= diagonal_temp;

        if(row_convergence(i, recvbuf, m_order, eps) == 0)
            loop = 0;
    }
   
    /*Process lines from the tail*/
    if(rank == source && tail != 0)
    {
        for(i = 0; i < lines_tail; i++)
        {
           
            index = m_order_pow_2 - (i*(m_order+1)) - 1;
            indexAdjust =  m_order_pow_2 - ((i+1)*m_order);
            diagonal_temp = a[index];
            a[index] = 0;
            int c = m_order - i - 1; 

            for(j = m_order - 1; j >= 0; j--)
            {
                if(j != c)
                    a[j+indexAdjust] /= diagonal_temp;
            }

            b[m_order - i - 1] /= diagonal_temp;

            if(row_convergence(m_order - i - 1, a, m_order, eps) == 0)
                loop = 0;
        }
    } 

    if(rank == source) {
        for(i = 0; i < m_order; i++)
        {
            x_res[i] = x[i];
        }
    }

    
    for(it = 0; it < max_it && loop == 1; it++)
    {
        /*Broadcast last Xk values to all tasks*/
        MPI_Bcast(x_res,m_order,MPI_DOUBLE,source,MPI_COMM_WORLD);
    
        for(i = 0; i < m_order; i++)
            xTmp[i] = x_res[i];

        for(i = 0; i < lines_per_task; i++)
        {
            index = (rank*lines_per_task) + i;
            indexAdjust = i*m_order;

            /*new*/
            dbuf[i] = xTmp[index];
            
            x_res[index] = bbuf[i];

            for(j = 0; j < m_order;j++)
            {
                   x_res[index] -= recvbuf[j+indexAdjust]*xTmp[j];
                   
            }
            
            dbuf[i] -= x_res[index];

       }

           /*Calculate X(k+1) for each line in the tail and stores*/ 
          if(rank == source && tail != 0)
          {
            for(i = 0; i < lines_tail; i++)
            {
                indexAdjust =  m_order_pow_2 - ((lines_tail - i)*m_order);
                int c = m_order - lines_tail + i;
                dif[c]  = xTmp[c];  
                x_res[c] = b[c];
                for(j = 0; j < m_order; j++)
                {
                        x_res[c] -= a[j+indexAdjust]*xTmp[j];
                }

                dif[c] -= x_res[c];
            }

        }

        //MPI_Allreduce(&max,&globalMax,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
        MPI_Gather( dbuf, lines_per_task, MPI_DOUBLE, 
            dif,lines_per_task, MPI_DOUBLE, source, MPI_COMM_WORLD);

        if(rank == source)
        {
           // printf("max %lf\n %d", max2(dif,x_res,m_order),it);
            max2(dif,x_res,m_order);
        }

       /* if(globalMax <= eps)
        {
            printf("It %d\n",it);
            loop = 0;
        }*/

        /*Save the result for each line in the root*/
        MPI_Gather( &x_res[rank*lines_per_task], lines_per_task, MPI_DOUBLE, 
            x_res,lines_per_task, MPI_DOUBLE, source, MPI_COMM_WORLD);


        MPI_Barrier(MPI_COMM_WORLD);


    }
     
        /*Get the values for A* and B* from other tasks*/
        MPI_Gather( recvbuf, recvcount, MPI_DOUBLE, 
            a,sendcount, MPI_DOUBLE, source, MPI_COMM_WORLD);
        MPI_Gather( bbuf, lines_per_task, MPI_DOUBLE, 
            b,lines_per_task, MPI_DOUBLE, source, MPI_COMM_WORLD);
        


        if(rank == source)
        {
            check_result(row_test,a,x_res,b, m_order,eps);
        }

    MPI_Finalize();
    return x;
}


double *solve(double *a, double *b, double *x, int m_order, int row_test,int max_it, double eps)
{

    unsigned int m_order_pow_2 = m_order * m_order;
    int temp = 0;
    int cond = 0;
    int j = 0;
    int i = 0;
    int it;
    int index;
    double *xTmp = malloc(m_order * sizeof(double));
    double *dif = malloc(m_order * sizeof(double));

    for(i = 0; i < m_order_pow_2; i+=(m_order+1))
    {
        temp = i - i%m_order;
        cond = temp + m_order;
        for(j = temp; j < cond; j++)
        {
            if(j != i)
            {
                a[j] /= a[i];
            }
        }
        b[i%m_order] /= a[i];
        a[i] = 0;
    }

    if(row_eval(a,m_order,eps) == 1)
        printf("Convergence test passed\n");
    else {
        printf("Convergence test failed\n");    
        return 0;
    }

    for(it = 0; it < max_it ; it++)
    {
        for(j = 0; j < m_order; j++)
        {
            dif[j] = xTmp[j] = x[j];
        }

        for(j = 0; j < m_order; j++)
        {
            x[j] = b[j];
            index = m_order*j; 
            for(i = 0; i < m_order; i++,index++)
            {
                x[j] -= (a[index]*xTmp[i]);
            }
            dif[j] -= x[j];
        }

        if(max(dif,x,m_order) <= eps)
            break;
    }
 
    check_result(row_test,a,x,b, m_order,eps);

    return x;
}

double max(double *vdif, double *vx, int size)
{
    double max,max2;
    int i;
    double temp,temp2;

    max = max2 = 0;

   for(i = 0; i < size; i++)
   {
        temp = fabs(vdif[i]);
        temp2 = fabs(vx[i]);
        if(temp > max)
            max = temp;
        if(temp2 > max2)
            max2 = temp2;
   }

   return max/max2;
}

int row_eval(double *matrix, int order, double eps)
{
    int i;
    for(i = 0; i < order; i++)
    {
        if(row_convergence(i,matrix,order,eps) == 0)
            return 0;
    }
    return 1;
}

static int row_convergence(int row, double *matrix, int order, double eps)
{
    int row_index = order * row;;
    int i = row_index;;
    double  acc = 0;
    double diff;

    do
    {
        if(i != (row_index + row))
        {
            acc += fabs(matrix[i]);
        }
        i++;
    } while(i < row_index + order);
    diff = 1 - acc;

    return diff > -eps;
}

static void check_result(int row, double *matrix, double *x, double *b, int order, double eps)
{
    int row_index = order * row;;
    int i = row_index;
    int j = 0;
    double  acc = 0;
    double diff;

    do
    {
        if(i != (row_index + row))
        {
            acc += matrix[i]*x[j];
        }else
            acc += x[j];
        i++;
        j++;
    } while(i < row_index + order);
    printf("%lf ?= %lf\n", acc,b[row]);

    //diff = fabs(b[row] - acc);
    //printf("Line %d tested = %lf\n",row,diff);
}
